```{r message=FALSE, warning=FALSE, paged.print=FALSE}
###########
# Setting #
###########
rm(list=ls(all=TRUE))
options("contrasts")
options(contrasts = c("contr.helmert", "contr.poly"))

library(tidyverse)
library(openxlsx)
library(car)
library(haven)
library(olsrr)
library(psych)
library(multcomp) 
library(scales)
library(gridExtra)
library(psych)

filter <- dplyr::filter
select <- dplyr::select
alpha <- psych::alpha

data <- read_sav("The Curve Not Taken_data.sav")
data <- as_tibble(as.data.frame(data)) %>% 
  filter(treatment != '')
```

```{r}
alpha(data[,c('HelAtti3', 'HelAtti4',
              'ExpertPer4', 'ExpertPer5',
              'polAtti2.9', 'polAtti2.8')])

alpha(data[,c('PolBias6','PolBias7', 'PolBias8')])
```

```{r}
#################
# Data cleaning #
#################
data$treatment <- as.factor(data$treatment)
data$factcheck <- as.factor(data$factcheck)
data$defense <- as.factor(data$defense)
data$racism <- as.factor(data$racism)
data$polParty <- as.factor(data$polParty)

mydata <- data %>% 
  mutate(healthExpertProfessionalism = 
           (PolBias6 + PolBias7 + PolBias8)/3,
         TestAndTrace3 = (HelAtti3 + HelAtti4 + 
                            ExpertPer4 + ExpertPer5 + 
                            polAtti2.8 + polAtti2.9)/6) %>%
  select(treatment, factcheck, 
         defense, racism, polParty, 
         Compa1, Compa3, Compa8,
         TestAndTrace3, Crespon2,
         healthExpertProfessionalism, 
         ResponseId)
```
       
              
```{r}
# Cook's distance function
cook.fuc <- function (model, data) {
cooksd <- cooks.distance(model)
threshold <- 4/(dim(data)[1])
outlier <- which(cooksd > threshold)
}
```

```{r}
#########################
# H1_a: Testing & Tracing
#########################

model11.data <- mydata %>% 
  drop_na(treatment, factcheck, 
          defense, racism, 
          polParty, Compa1)

model11 <- lm(Compa1 ~ treatment + 
                factcheck + defense + 
                racism + polParty +  
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty,
              data = model11.data)

model11.outlier <- cook.fuc(model11, model11.data)
model12.data <- model11.data[-model11.outlier,]

model12 <- lm(Compa1 ~ treatment + 
                factcheck + defense + 
                racism + polParty +  
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty,
              data = model12.data)

model12.result <- Anova(model12, type = 3)
model12.result
```


```{r}
# Tukey-HSD for H1_a

model12.aov <- aov(Compa1 ~ treatment + 
                factcheck + defense + 
                racism + polParty +  
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty,
              data = model12.data)

TukeyHSD(model12.aov)$`treatment:polParty`
```

```{r}
##############################
# H1_b: Lockdown misperception
##############################

model21.data <- mydata %>% 
  drop_na(treatment, factcheck, 
          defense, racism, 
          polParty, Compa3)

model21 <- lm(Compa3 ~ treatment + 
                factcheck + defense + 
                racism + polParty +  
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty,
              data = model21.data)

# Cook's distance for removing outliers
model21.outlier <- cook.fuc(model21, model21.data)
model22.data <- model21.data[-model21.outlier,]

model22 <- lm(Compa3 ~ treatment + 
                factcheck + defense + 
                racism + polParty + 
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty,
              data = model22.data)

model22.result <- Anova(model22, type = 3)
model22.result
```

```{r}
# Tukey-HSD for H1_b

model22.aov <- aov(Compa3 ~ treatment + 
                factcheck + defense + 
                racism + polParty + 
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty,
              data = model22.data)

TukeyHSD(model22.aov)$`treatment:polParty`
```

```{r}
#######################################
# H1_c: Policy similarity misperception
#######################################

model31.data <- mydata %>% 
  drop_na(treatment, factcheck, 
          defense, racism, 
          polParty, Compa8)

model31 <- lm(Compa8 ~ treatment + 
                factcheck + defense + 
                racism + polParty +  + 
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty, 
              data = model31.data)

# Cook's distance for removing outliers
model31.outlier <- cook.fuc(model31, model31.data)
model32.data <- model31.data[-model31.outlier,]

model32 <- lm(Compa8 ~ treatment + 
                factcheck + defense + 
                racism + polParty + 
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty,
              data = model32.data)

model32.result <- Anova(model32, type = 3)
model32.result
```

```{r}
####################
# H2: Policy support
####################

model41.data <- mydata %>% 
  drop_na(treatment, factcheck, 
          defense, racism, 
          polParty, TestAndTrace3)

model41 <- lm(TestAndTrace3 ~ treatment + 
                factcheck + defense + 
                racism + polParty + 
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty,
              data = model41.data)

# Cook's distance for removing outliers
model41.outlier <- cook.fuc(model41, model41.data)
model42.data <- model41.data[-model41.outlier,]

model42 <- lm(TestAndTrace3 ~ treatment + 
                factcheck + defense + 
                racism + polParty +  
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty, 
              data = model42.data)

model42.result <- Anova(model42, type = 3)
model42.result
```

```{r}
########################
# H3: Presidential Blame
########################

model51.data <- mydata %>% 
  drop_na(treatment, factcheck, 
          defense, racism, 
          polParty, Crespon2)

model51 <- lm(Crespon2 ~ treatment + 
                factcheck + defense + 
                racism + polParty + 
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty, 
              data = model51.data)

# Cook's distance for removing outliers
model51.outlier <- cook.fuc(model51, model51.data)
model52.data <- model51.data[-model51.outlier,]

model52 <- lm(Crespon2 ~ treatment + 
                factcheck + defense + 
                racism + polParty + 
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty, 
              data = model52.data)

model52.result <- Anova(model52, type = 3)
model52.result
```

```{r}
#############################
# H4: Trust in health experts
#############################

model61.data <- mydata %>% 
  drop_na(treatment, factcheck, 
          defense, racism, 
          polParty, healthExpertProfessionalism)

model61 <- lm(healthExpertProfessionalism ~ treatment + 
                factcheck + defense + 
                racism + polParty + 
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty,
                    data = model61.data)

# Cook's distance for removing outliers
model61.outlier <- cook.fuc(model61, model61.data)
model62.data <- model61.data[-model61.outlier,]

model62 <- lm(healthExpertProfessionalism ~ treatment + 
                factcheck + defense + 
                racism + polParty +  
                treatment * factcheck + 
                treatment * defense + 
                treatment * racism + 
                treatment * polParty, 
                    data = model62.data)

model62.result <- Anova(model62, type = 3)
model62.result
```

```{r}
#########
# Graph #
#########
myrecord.treat <- c('Control', 'Outcome', 'Policy')
myrecode.party <- c('Republicans', 'Democrats', 'Independents')

# H1_a
graph1 <- model12.data %>%
  select(treatment, polParty, Compa1) %>%
  group_by(treatment, polParty) %>%
  summarise(Compa1 = mean(Compa1, na.rm =T))

myplot1 <- ggplot(graph1, aes(x=myrecord.treat[treatment], 
                              y=Compa1, 
                              group = myrecode.party[polParty])) + 
  geom_line(aes(linetype=myrecode.party[polParty]))+ 
  geom_point(aes(shape=myrecode.party[polParty])) + 
  scale_linetype_manual(values=c("solid","twodash", "dotted"))

mydisplay <- theme(panel.background = element_rect(fill='grey90',
                                                   color = NA),
                   axis.title.x = element_blank(),
                   legend.title = element_blank(),
                   axis.title.y = element_text(face="bold"),
                   text=element_text(family="Arial"))

lab1 <- labs(y = "Testing & Tracing")
deci <- scale_y_continuous(labels = number_format(accuracy = 0.1, 
                                    decimal.mark = '.'))

M1.plot <- myplot1 + mydisplay + lab1 + deci
M1.plot
```

```{r}
# H1_b
graph2 <- model22.data %>%
  select(treatment, polParty, Compa3) %>%
  group_by(treatment, polParty) %>%
  summarise(Compa3 = mean(Compa3, na.rm =T))

myplot2 <- ggplot(graph2, aes(x=myrecord.treat[treatment], 
                              y=Compa3, 
                              group = myrecode.party[polParty])) + 
  geom_line(aes(linetype=myrecode.party[polParty]))+ 
  geom_point(aes(shape=myrecode.party[polParty])) + 
  scale_linetype_manual(values=c("solid","twodash", "dotted"))

lab2 <- labs(y = "Lockdown Misperception")

M2.plot <- myplot2 + mydisplay + lab2 + deci
M2.plot
```

```{r}
# H1_c
graph3 <- model32.data %>%
  select(treatment, polParty, Compa8) %>%
  group_by(treatment, polParty) %>%
  summarise(Compa8 = mean(Compa8, na.rm =T))

myplot3 <- ggplot(graph3, aes(x=myrecord.treat[treatment], 
                              y=Compa8, 
                              group = myrecode.party[polParty])) + 
  geom_line(aes(linetype=myrecode.party[polParty]))+ 
  geom_point(aes(shape=myrecode.party[polParty])) + 
  scale_linetype_manual(values=c("solid","twodash", "dotted"))

lab3 <- labs(y = "Policy Similarity Misperception")

M3.plot <- myplot3 + mydisplay + lab3 + deci
M3.plot
```

```{r}
# H2
graph4 <- model42.data %>%
  select(treatment, polParty, TestAndTrace3) %>%
  group_by(treatment, polParty) %>%
  summarise(TestAndTrace = mean(TestAndTrace3, na.rm =T))

myplot4 <- ggplot(graph4, aes(x=myrecord.treat[treatment], 
                              y=TestAndTrace, 
                              group = myrecode.party[polParty])) + 
  geom_line(aes(linetype=myrecode.party[polParty]))+ 
  geom_point(aes(shape=myrecode.party[polParty])) + 
  scale_linetype_manual(values=c("solid","twodash", "dotted"))

lab4 <- labs(y = "Policy Support")

M4.plot <- myplot4 + mydisplay + lab4 + deci
M4.plot
```

```{r}
# H3
graph5 <- model52.data %>%
  select(treatment, polParty, Crespon2) %>%
  group_by(treatment, polParty) %>%
  summarise(Crespon2 = mean(Crespon2, na.rm =T))

myplot5 <- ggplot(graph5, aes(x=myrecord.treat[treatment], 
                              y=Crespon2, 
                              group = myrecode.party[polParty])) + 
  geom_line(aes(linetype=myrecode.party[polParty]))+ 
  geom_point(aes(shape=myrecode.party[polParty])) + 
  scale_linetype_manual(values=c("solid","twodash", "dotted"))

lab5 <- labs(y = "Presidential Blame")

M5.plot <- myplot5 + mydisplay + lab5 + deci
M5.plot
```

```{r}
# H4
graph6 <- model62.data %>%
  select(treatment, polParty, 
         healthExpertProfessionalism) %>%
  group_by(treatment, polParty) %>%
  summarise(healthExpertProfessionalism = 
              mean(healthExpertProfessionalism, na.rm =T))

myplot6 <- ggplot(graph6, aes(x=myrecord.treat[treatment], 
                              y=healthExpertProfessionalism, 
                              group = myrecode.party[polParty])) + 
  geom_line(aes(linetype=myrecode.party[polParty]))+ 
  geom_point(aes(shape=myrecode.party[polParty])) + 
  scale_linetype_manual(values=c("solid","twodash", "dotted"))

lab6 <- labs(y = "Trust in Health Experts")

M6.plot <- myplot6 + mydisplay + lab6 + deci
M6.plot
```
